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Cluster Monte Carlo Algorithms for Dissipative Quantum Systems 

Philipp Werner and Matthias Troyer 
Institut fiir theoretische Physik, ETH Zurich, CH-8093 Zurich, Switzerland 

We review eflicient Monte Carlo metiiods for simulating quantum systems which couple 
to a dissipative environment. A brief introduction of the Caldeira-Leggett model and the 
Monte Carlo method will be followed by a detailed discussion of cluster algorithms and 
the treatment of long-range interactions. Dissipative quantum spins and resistively shunted 
Josephson junctions will be considered. 

§1. Introduction 

Dissipation and decoiierence in quantum systems liave important practical im- 
plications, from stabilizing superconductivity in granular materials, to the loss of 
information stored in qubits. In all the models we study, energy dissipation is in- 
troduced by coupling a quantum mechanical degree of freedom to a bath of suitably 
chosen harmonic oscillators. The degrees of freedom of this bath can be integrated 
out analytically, which leads to an effective action with long-range interactions in 
imaginary time. In this article we will review and present a number of recently 
developed efficient cluster Monte Carlo algorithms for various dissipative quantum 
systems. These algorithms have been used to study spin chains with dissipation 
coupling to the site variables and resistively shunted Josephson junctions. ^^^'^^ We 
will also present relevant theoretical background for these models and mention the 
open questions which have been addressed in the numerical investigations. 

§2. Dissipative quantum models 

2.1. Caldeira-Leggett model 

While the concept of coupling a quantum mechanical degree of freedom to an 
environment goes back to the early days of quantum mechanics, it was the work of 
Caldeira and Leggett^^'^^ which pointed out the important consequences of dissipa- 
tion at a time when experimental results on macroscopic quantum tunneling became 
available. To introduce their model we will consider a classical system whose dissi- 
pative dynamics in terms of some coordinate q is described in the absence of external 
forces by the phenomenological equation of motion 

Mq{t)+m{t) + [V'{qm=^, (2-1) 

where ry denotes the friction coefficient. Because the frictional force is proportional 
to q, one calls the dissipation "Ohmic". 

The equation of motion (|2-ip can be obtained from a microscopic model by 
coupling the system linearly in the coordinates to an appropriately chosen set of 
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harmonic oscillators, that is, from the Lagrangian 
M 



L 



+ \^\ ^c,il - rriaujl (xa ^9 ) \ , (2-2) 



where Xq, ma, oJa and Cq denote the oscillator positions, masses, frequencies and 
coupling strengths, respectively. Defining the spectral density 

J{u:):=^-Y,^5{u:-Ua), (2-3) 
2 ^-^ rriniVn, 

we can reproduce the phenomenological equation (|2-1|) by setting 

J{uj) = r]u;, (2-4) 

for oj > O.^^'^)'^''-* In any real physical system, the harmonic oscillators representing 
the environment will have finite frequencies, so the linear relationship (|2-4|) must be 
cut off at some value ujq. 

The partition function of the quantum mechanical system can be obtained from 
the classical Euclidean Lagrangian 



P M n 



ruaxl + niaiJi ( Xa ) > , (2-5) 



by the path integral formalism.^)' ^^'^'^^ Integrating out the harmonic oscillator de- 
grees of freedom leads to an effective action 



/\\2 



Jo (sin(||r-T'|))2- 

(2-6) 

In analytical calculations, but also for certain algorithms discussed in chapter 3, 
it is useful to express the damping term (proportional to ij) in Fourier space. In the 
limit P ^ oo and for asymmetric Fourier transforms one obtains 



oo 



'S'damp = ^ y duj\uj\\q{uj)\'' , (2-7) 

which is local in the frequency lo, in contrast to the imaginary-time integral p-6j) . 

2.2. Dissipative Quantum Spin Chains 

The effective action for a single 0{n) spin with on-site dissipation is obtained 
from Eq. (|2-6() by choosing an n-component vector q as the coordinate which couples 
to the environment and a double-well potential 

V{q) = '-q' + ^/ (2-8) 
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with minima at \q\ = 1. Coupling the spins spatially through nearest neighbor bonds 
of the form —Kijj ■ qj-\-i then yields the action for the spin chain 



dT 







47r ^ 



— {drqj) - Kqj ■ Qj+i + -qj + —qj 



:r,r,r'(^y^Mp^. ,2.9) 
„J„ \fll (sin(J|T-T'|))2 ^ ' 



The special cases n = 1 (Ising spin chain) and n = 2 (XY-spin chain) will be discussed 
in chapter 4. 

To simplify the calculations, we will discretize the action ()2-9|) in imaginary time, 
noting that the cutoff of the linear spectral density J{uj) at uJc would anyhow lead to 
modifications in /Sdamp at small t — t' . We use a Trotter number A',-, corresponding 
to a discretization step At = fi/Nr- Furthermore we employ a "hard spin" approx- 
imation and choose the potential ()2-8() such that the spin states at any time can 
be described as 0{n) vectors a of norm \a\ = 1. Introducing the the coefficients 
r = Mj At and a = r//vr, the discretize action becomes 

^ = - E E • ^^+1-^ + ^^^^^ ■ - "E E (sinS-V-I'i))2-' 

j=i fc=i j=ik<k' ^ ^ ^^^ri 

(2-10) 

which can be simulated as a two-dimensional system of classical 0(n)-spins with 
short-range interactions in the "space"- and long-range interactions in the "imaginary 
time" direction. 

2.3. Resistively shunted Josephson junction 

Another example of a dissipation coupled quantum system is a Josephson junc- 
tion with Josephson coupling energy Ej and capacitance C, which is shunted by an 
Ohmic resistor Rg , as illustrated in Fig. ^ While Cooper pairs can tunnel through 
the junction, which consists of a thin layer of insulating material or a constriction in 
the superconductor, electrons can flow through the resistor and dissipate energy. Let 
and (j)r denote the phases of the superconducting order parameter on the left and 
right island, respectively, and (p = (pi — (pr the phase difference across the junction. 
This phase difference is related to the voltage drop V across the junction and the 
super-current Is by the Josephson relations {h/2e is the flux quantum) 

y = A^, (2.11) 

Is = lcsm4>, (2-12) 
where Ic denotes the critical current, which is related to the coupling energy by 

Ej = ^. (2-13) 
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Fig. 1. Illustration of a resistively shunted Josephson junction. It consists of two superconducting 
islands connected by a Josephson junction of coupling energy Ej and capacitance C and a 
parallel resistor Rs through which electrons can flow 



The dynamics of such a device can be determined from a balance of currents, which 
for zero external current reads 

V 

— -^displacement ~l~ -^shunt ~l~ -^supra — CV H — h 1^ sin (j) 



hc^ ^ I hc^ 2eEj . 

= C — sm 6. 

2e Rs2e h ^ 

If we identify 5^ with a coordinate q and introduce the potential 



Uiq) = -Ejcosi(l)iq)) = - E j cos{{2e / fi)q) , 



then ()2-14|) becomes 



Cq+^q + U'{q) 
Us 



0, 



(2-14) 



(2-15) 



(2-16) 



which is just equation (|2-1|) with the substitutions M ^ C and V ^ According 
to the recipe of Caldeira and Leggett outlined in section 2.1, the dissipation in the 
shunt resistor can be reproduced by coupling q to an Ohmic heat bath, which is 
integrated out and leads to the effective action (|2-6() . Expressing g, M and rj again 
in terms of the original variables C and Rs, one finds the action for the resistively 
shunted Josephson junction. Introducing furthermore the effective charging energy 

2 

of the junction Ec = ^ (which sets the energy scale) and the quantum of resistance 
Rq = -Aj, the action becomes (setting h = 1) 



S[<p] 



13 



dT 



1 



1 



16Ec Vdr/ 
Rq 



E J cos( 



87r2 R, 



JO 



drdr' 



sin((7r//3)(r-T'))2 



(2-17) 



Resistively shunted Josephson junctions at zero temperature undergo a super- 
conductor-to-metal transition if the shunt resistance Rg is increased beyond the 
critical value Rq. This dissipative phase transition was first predicted by Schmid^^-* 
and Bulgadaev^^^ and subsequently studied by several authors. ^^^"^''^ 
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Fig. 2. Circuit diagram of the single electron box. The box with excess charge n is indicated by 
the dashed rectangle. It is connected to a voltage source through a capacitor Co and a tunnel 
junction with resistance Rt and capacitance Ct- 



2.4. Single electron box 

The so-called single electron box, which consists of a low-capacitance metallic 
island connected to an outside lead by a tunnel junction, is described by the action 
of a dissipative quantum rotor. Due to the large charging energy of the island, the 
presence of excess charges influences single electron tunneling and such a device thus 
exhibits Coulomb blockade phenomena. 

A circuit diagram of a single electron box is shown in Fig. [2 The box with excess 
charge n is controlled by an external voltage source Vg, to which it is connected 
through a capacitor Cq and a tunnel junction with resistance Rt and capacitance 

2 

Ct- The bare charging energy Ec = 2{Ct+CG) ^^^^ energy scale. An applied gate 
voltage Vg induces a continuous polarization charge ug = CcVc/e. 

The effective action of the single electron box with gate charge zero can be 
derived from a microscopic theory, as detailed in Refs. 7), 16), 18)-20): 

" 4EcJo 2m^Rt Jo (siii((ir//3)(T - t')))^ ^ ' 

The compact angular variable 9 is conjugate to the number of excess charges on 
the island and the partition function can be written as the sum over all paths with 
winding number lo = in, n = 0, 1, 2 . . .. 

At low temperatures, ksT <^ Ec, and high tunneling resistance, Rt » Rr = 
hje^ ~ 25.8A:i7, the number of excess charges on the island is a staircase function 
with unit jumps at ug = 1/2 (mod 1). Thermal fluctuations will round off the corners 
of the step structure. But even at zero temperature, electron tunneling processes 
can smear out the staircase, which for large tunneling approaches a straight line. 
These quantum fluctuations renormalize the ground state energy and lead to an 
effective charging energy E^^, which can be computed from the expectation value of 
the winding number oj squared, 

E*c _ 27r2 , 27r2 1 f^^ ^ , /■e„±2.n 



Ec PEc 



pEc ^ Jo Jdo 



In chapter 4 we will address the behavior of E^/Ec in the limit of large tunneling 
conductance. 
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§3. Algorithms 

In this main section we present algorithms used for the simulation of the dissi- 
pation coupled quantum systems discussed in the previous chapter. A brief intro- 
duction to the Monte Carlo method and the concept of importance sampling will 
be followed by a section on cluster updates for spin systems and the efficient treat- 
ment of long-range interactions. Cluster algorithms identify clusters of spins which 
can be updated simultaneously in a rejection-free move. Near phase transitions, 
where domains of aligned spins occur on all length scales, they dramatically reduce 
autocorrelation times (the number of updates needed to produce an uncorrelated 
configuration) compared to local update schemes. 

3.1. Importance sampling and the Metropolis Algorithm 

We would like to compute the expectation value of some observable O, 



(0) = j;P(x)0(x), (3-1) 



n 

where P denotes a probability distribution defined on the configuration space Q. 
Since it is usually not possible to sum over all configurations, the strategy is to 
sample random configurations xi, . . . , x^r, according to their importance P{xi) and 
approximate Eq. (|3-1|1 by 

_ 1 ^ 

^^]^E0(^*)- (3-2) 

■t=i 

For statistically independent 0(xj) and large enough N it follows from the central 
limit theorem that O follows a Gaussian distribution centered at (O) and standard 
deviation 

_ /Var(O) ja^-O'^ , , 

As a consequence, the accuracy of a Monte Carlo simulation scales as Xj^ff^ with 
computation time N . 

The challenge is to generate a chain of configurations xi, . . . , x^r such that the 
Xi are sampled according to their weight P{xi). For this purpose we introduce a 
Markov process which defines a probability distribution Pj+i at time (Monte Carlo 
step) t + 1 from the distribution Pt at time t by means of a transition matrix T, 

Pt^r{y) = Y,Ty\M^). (3-4) 

X 

The elements Ty\x of T correspond to the transition probabilities from state x to 
state y and therefore 

>0, J^^Tyi, = 1. (3-5) 

y 

As a consequence -Pt+i is a probability distribution if Pf is one. An additional 
requirement for T is ergodicity: it must be possible to reach every configuration y 
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from any configuration a; in a finite number of Markov steps. This assures that in 
the limit of infinitely many steps, the whole configuration space is sampled. Finally, 
T must be such that the probability distribution P is stationary, 

P{y) = Y,TyixP{x). (3-6) 

X 

This condition is satisfied, if 

Ty\,P{x) = T,\yPiy), (3-7) 

as is obvious by summing over x. The relation 1)3 Tj) is referred to as detailed balance. 

Given the above, any probability distribution Po{x) = 6x,xi will converge to P{x) 
as t ^ oo.^^) Hence, an infinitely long Markov chain starting with Xi and defined by 
the transition matrix T generates configurations according to their weight P{x). 

We would now like to explicitly construct the transition matrix T. To this 
end we decompose the transition probability to move from a configuration x to a 
configuration y, Ty^^, into an 'a priori' probability P^^^ to propose the move from 
X to y and pf?!; to accept the move. Assuming that p^T"^ = and satisfies 

y\x x\x y\x 

conditions ()3-5() . the transition matrix defined as 

satisfies Eq. (|3-5|) as well. The detailed balance condition 1)3 -71) now reads 

p7J'pW^nx)=p':\^pT^p{y), (3-9) 

and can be satisfied in several ways. The Metropolis algorithm^^) is based on the 
choice 

P^|.^ = --(^l'^Spp(^J- (3-10) 

3.2. Cluster algorithms for classical spins 

The Monte Carlo simulation of spin systems near criticality requires the use 
of efficient algorithms, which we will first discuss for the Ising model with nearest 
neighbor interactions {K > 0, denotes a pair of nearest neighbor sites) 

S = -K^aiaj. (3-11) 

Since the correlation length diverges in the vicinity of a phase transition, large do- 
mains of aligned spins appear. Thus, local update schemes become very inefficient 
in changing the overall configuration, which leads to large autocorrelation times. In 
order to overcome this problem, which means to induce transitions between probable 
configurations which differ on large scales, it is necessary to fiip entire clusters of 
spins in each Monte Carlo update. The Swendsen-Wang algorithm^^^ chooses these 
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spin rep. C„ C„+i 

joint rep. (C„,G„) ---> (C„+i,G„) ---> (C„+i,G„+i) 

\ / rG„+i|c„+i 
bond rep. Gn 

Fig. 3. Illustration of the Markov process in the Swendsen-Wang cluster algorithm, which switches 
back and forth between spin configurations (Cn) and bond configurations (G„). The transition 
probabilities rc„+i|G„ and rG„+i|c„+i can be found from the joint representation in terms of 
spin and bond variables. 



clusters in such a way, that detailed balance is automatically satisfied (no cluster 
move is rejected). 

The Swendsen-Wang algorithm is an example of a dual Monte Carlo algo- 
rithm,^^) which switches back and forth between two configuration spaces, as illus- 
trated in Fig. |3J One representation is in terms of Ising spin variables ai G {—1, 1}, 
the other one is the random cluster representation of Fortuin and Kasteleyn,^^) which 
is expressed in terms of bond variables bij S {0, 1}. The transition probabilities be- 
tween configurations in both spaces follow directly from the joint Edwards-Sokal 
representation in terms of spin and bond variables. 

The partition function of the Ising model ()3-ll() in these three representations 
reads 

• spin representation 

Z = Y^ e^^<».^> '''''' =Y.]1 e^''^''^ (3-12) 

• joint representation 

^ = E E n Koe-'^ + ^^...1^'^.'^. (1 - ^-"'')] ' (3-13) 

{b} W} (ij) 



bond representation 



^^g-K7vr^(e2i^_l)n,2nc^ (3.14) 

{b} 



where Nj^ denotes the total number of bonds, Uh the number of bond variables 



with value bij = 1 and Ur. the number of clusters. 



The joint representation in Eq. 1)3- is obtained from Eq. H3-12() by introducing 
bond variables bij £ {0, 1} using the identity 



= E K-oe-^^ + - . (3-15) 



Summing over the spin variables in Eq. ()3-13() then yields the Fortuin-Kasteleyn 
random cluster representation, Eq. H3-14() . 
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The Swendsen-Wang algorithm follows directly from the joint representation.^'''' 
Given a spin configuration {crj}, the conditional probability for freezing a bond 
{bij = 1) is found from Eq. H3-13() to be 



1 - p-'^^ 

^ 6 - -2K 



Tfe^,=i|.,^., = 0. (3-17) 
If the bond configuration {bij} is given, we find 

Tcj,=aj\h,=i = 1, (3-18) 
r.,=.,|6,,=o = l/2, (3-19) 

which means that all spins of a given cluster point in the same direction, whereas 
the spins on different clusters are uncorrelated after the update. 

In the more general case, where spins interact over arbitrary distances with 
different strengths, the probability of a bond between the sites i and j becomes 

nond(^,j) =max(0,l-exp[-Z\S(i,i)]), (3-20) 

where AS{i,j) = 5'updated(^j j) — 'S'originai(*, j) denotes the cost in action of fiipping a 
spin associated with the pair 

The Swendsen-Wang algorithm proceeds in the following three steps, which for 
the case of the dissipative quantum Ising chain ()2-10|) are illustrated in Fig. 0] 

1. Bonds are inserted between the sites of interacting spins with probability ()3-2U|) . 

2. Clusters of connected spins are identified. 

3. The spins of each cluster are fiipped with probability ^. 

Wolff ^^-^ extended the ideas of Swendsen and Wang to 0{n) spin models. First of 
all he noted that it is more efficient to build a single cluster from a randomly chosen 
site and update all the spins in this cluster. This has the advantage that most 
updates take place in the largest regions of aligned spins. By choosing a random 
direction e on the n-sphere and projecting the spins onto this direction, one ends 
up with an Ising like system. Updating a spin then means fiipping this projected 
component, or in other words mirroring the spin on the plane perpendicular to e 

a ^ a - 2{a ■ e)e. (3-21) 

The probability (|3-2U|) for inserting a bond between two spins at sites i and j remains 
valid. A Wolff cluster update therefore proceeds in the following three steps 

1. Choose a random site and a direction e. 

2. Find all spins connected to this site using the bond probability (|3-2(jp . where 
AS{i,j) = S{ai- 2((Ti • e)e,aj) - S{ai,aj) 

3. Mirror all the spins in the cluster on the plane perpendicular to e. 

3.3. Efficient treatment of long-range interactions 

In the presence of long-range forces, every site interacts with every other site 
and an algorithm which iterates over all bonds would be of order 0(A^^), where 
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Fig. 4. Illustration of the Swendsen-Wang cluster algorithm. The figures show clockwise from the 
top left: (i) original spin configuration, (ii) insertion of bonds, (iii) clusters of connected spins, 
(iv) new spin configuration obtained by flipping the spins of each cluster with probability ^. 
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N is the number sites. This is prohibitively slow and thus we present here two 
more sophisticated algorithms, proposed by Luijten and Blote,^*^-* which are of order 
O(A^logA^) and 0{N) respectively. These ideas are essential for the simulation of 
dissipative quantum systems at low temperature, since we have seen in chapter 2 
that dissipation introduces long-range interactions in the imaginary time direction. 
We consider first a chain of classical Ising spins with long-range interactions, 

5 = j>,aj. (3-22) 
The kernel g of interest in our applications is of the form (j ^ 0) 



(Mm' 



5(i)=«T^^7^^> (3-23) 



but this is of no importance for the following discussion. Prom Eq. (|3-2nj) it follows 
that the probability of a long-range bond between two parallel spins at sites and j 
is 

Pbond(0;j) = l-e-2^(^). (3-24) 

Hence, the probability that no bond is formed between spin and spins j'+l, . . . , n— 1 
(all assumed parallel) becomes 

n— 1 r n— 1 

^'nobond(0;j + l,...,n-l) = H e-29»=exp -2 ^ g{i) . (3-25) 

i=j+l I i=j+l 

If we define an array ("lookup-table") A of length A'^ with elements ^[0] = 1, and 
A[n] = nr=i e~'^9{i) forn = 1, . . . , iV - 1, then Eq. can be written as 

Pno bond (0;i + l,..., n-1) = ^^"""Z^ . (3-26) 

The spins connected to site are then calculated as follows: a random number 
ri is chosen in the interval [0, 1) and the array A searched for the first index, say 
rei, such that A[ni] < ri. For the calculation of the next bond, the values of A 
must be divided by A[ni], or equivalently, the random number r2 multiplied by 
^[ni]. The site number 712 of the second connected spin is the first index, such that 
A[n2] < j4[ni]r2. This process is continued until ^[n^.ijrfc < A[N — 1]. 

The outlined procedure assumes that all spins are aligned. Before actually in- 
serting a proposed bond, we have to check whether the two spins are parallel. If 
the bisection method is used to search for the lowest index such that A[n] < r, the 
calculation time is of order O(logA^) per site. An update of the whole lattice with 
lookup-table is therefore of order 0(A^log A^). 

The search of a lookup-table could even be avoided, if it is possible to evaluate 
Eq. H3-25|) analytically, as is the case for the kernel H3-23(l . To this end we approximate 
the sum in the exponent by an integral 

El 2 {j^Y vr / vr / 1\\ / / . 1 

KiK, a \ dT- r-L — --77 = —ol — cot — n — cot — ? H — 
y,+ i (sin(fr))2 iv[ \N\ 2)) KnV 2 



(3-27) 
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Given the value of j (index of the previous bond), the value of n (index of the 
next bond) is calculated from a uniformly distributed random number r E [0, 1) by 
inverting 



exp 



n-l 

Y. 

i=j+i 



1 — r, 



(3-28) 



where it is understood that n will be truncated to the next smaller integer (because 
the sum is only up to n — 1). Replacing the sum in Eq. ()3-28() by Eq. ()3-27|) we get 



n = int 



IN 1 

— I atan2 1, ^ h 

2 vr I ' tan(f (j + \)) vr 



iVln(l 



2a 



(3-29) 



The function atan2 satisfies tan(atan2(2;, y)) = x/y. It is defined on the interval 
[— vr, vr] and the signs of the two arguments are used to determine the quadrant. 

If the analytical formula (|3-29() is used to calculate bonds, the computation time 
per site is of order 0(1). In fact, from Eqs. (|3-25() and H3-27() it follows with j = 
and n = N, 



Pno bond = exp <^ 2a 



vr 

N 



cot ( vr I — cot ( I 

2NJ \2NJ 



-2a 



(3-30) 



where in the last step we assumed large. Hence, irrespective of the value of N, 
the probability for the algorithm to terminate at a given step is at least e~^" and 
thus the expected number of bonds is finite. 

Since short-range interactions are not well approximated by Eq. (|3-27j) . non- 
universal quantities, such as critical couplings, will have different values in this an- 
alytical implementation. However, universal quantities, such as critical exponents, 
will be calculated correctly, as they only depend on the asymptotic behavior of the 
interaction, which is accurately reproduced by Eq. H3-27() . 

Although the algorithm with analytical determination of bonded sites scales as 
0{N) and the algorithm based on a lookup table only as O(A^logA^), it turns out 
that the latter is faster by a factor of 2 for accessible system sizes. 

In the case of 0(n)-spins, where the bond strength depends on the spin compo- 
nents projected on some random direction e. 



rProj 



a ■ e, 



(3-31) 



one proposes bonds as explained for the Ising case. To account for the fact that the 
interaction is weaker than assumed in the look-up table, a random number r G [0, 1) 
is chosen and a proposed bond between sites i and j is only inserted if 



r(l 



< 1 



-29(«-i)o-. 



proj proj 



(3-32) 



In other words, bonds between 0(n)-spins are inserted with probability P]!^}"]^{h j) 
by proposing bonds with the higher probability -fbon|(^)i) and accepting them with 
probability P,^,^ (z, (z, j) . 
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3.4. Cluster algorithm for resistively shunted Josephson junctions 

In section 2.3 we derived the effective action for a resistively sliunted Josephson 
junction, which in terms of the dimensionless parameter a = Rq/Rs reads 



r-/3 

S[cP] = I dr 







l6Ec \dT 



Sjcos(0) 



in^Jo Jo sin((7r//3)(T-T'))2 ' 

(3-33) 

Because of the Josephson couphng energy and the non-compact nature of the phase 
variable in this action, one cannot directly employ the cluster algorithms available 
for spin systems, which were discussed in the previous section. In order to simulate 
resistively shunted junctions efficiently, we developed a new algorithm,^) consisting 
of two kinds of updates: (i) local updates in Fourier space compatible with the 
Gaussian terms in Eq. H3-33p and (ii) rejection-free cluster updates. The first type of 
moves assures ergodicity of the algorithm and the second type produces global cluster 
updates compatible with the energetic constraints from the Josephson potential. 

3.4.1. Local updates in Fourier space 

For the Monte Carlo simulation, we discretize imaginary time into N (assumed 
odd) time steps. The action (|3-33() can then be expressed in the simple form 

N-l N-1 

S[4>] = «A:|4l' - EjAt cos(0,-), (3-34) 

k=0 j=0 

where (f>k = Ylf=o e'^~^^4'j denotes the Fourier transform of (f). The positive coeffi- 
cients ak are defined as at = jj{go — Slk), with the Fourier transform of the kernel 

(i#o) 

Since cpl = 4>N-k-: only {0fc|A; = 0, . . . , {N — l)/2} need to be considered. 

In a local update of the frequency components (p^ and (pN-k: we choose a new 
value according to the probability distribution of the Gaussian term in Eq. (|3-34j) , 

~e-2»'=l'^fel', (3-36) 

by fixing the phase at a random value in the interval [0, 27r] and using exponen- 
tially distributed random numbers with mean l/(2afc). This move is accepted with 
probability 

p{<Pold ^ 0new) = min f 1, e-^^-^[<^— ^ (3.37) 



where iS'j[0] = —EjAT'^^Sf^cos{cl)j) and ipoid, 4'new denote the backward Fourier 
transform of the old and new /c-space configuration, respectively. Such local updates 
can be performed in a time 0{N) and have recently been used in the simulation of 
2D Josephson junction arrays. ^^-^ 
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Fig. 5. Illustration of local updates in Fourier space. Phase configurations = 0, . . . , A'^ — 1} 

are plotted as small circles and the value of the Josephson potential is sketched on the left. 
Dashed lines indicate the positions where the potential is minimal. The figures show from top 
to bottom: (i) Original phase configuration. The phase variables remain most of the time near 
a minimum of the cosine-potential, which results in a step-like structure, (ii) Randomly picked 
Fourier component, for which a new amplitude will be chosen according to Eq. H3-36|l . (iii) New 
phase configuration obtained if the new amplitude is close to zero. Since many of the phase 
variables end up in the region of high potential energy, this move will be rejected with high 
probability. 
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3.4.2. Cluster updates 

For reasonably large values of Ej, local fc-space updates which introduce phase 
changes on the order of 27r will be strongly suppressed, because their sinusoidal shape 
does not resemble an optimal phase slip path (see illustration in Fig. Algorithms 
based on local updates alone will therefore be ineffective near the phase transition, 
where phase slips start to proliferate. A typical path will stay most of the time 
near one or the other of the minima in the cosine-potential, as shown in Fig. El 
A simple idea for a global update compatible with this overall structure would be 
step- updates which shift the phases by ±27r in some random interval [jmim jmax] • 
It would, however, be better to let the algorithm choose itself the phase variables 

which can be shifted to a different potential minimum. 

The observation which leads to such a cluster algorithm is that a shift by mul- 
tiples of ib27r is not the only operation which leaves Sj invariant. The same is true 
for reflections on {mr}n£Z, that is, the positions of the maxima and minima of the 
cosine potential. We exploit the latter symmetry to design a rejection-free cluster 
update consisting of the following four steps, which are illustrated in Fig. |SJ 

1. An axis (p = '^'""^Tr with integral n^^^^ in the range [— nmax; '^max] is randomly 
chosen (the significance of nmax is discussed below) and a random site jroot is 
picked as the root site of the cluster. We introduce relative coordinates 

^'^^^ = - n'^'^vr, (3-38) 

which are updated in a cluster move as (j)^^^^ —(pj^^^, in complete analogy to 
the projected spin components H3-31() in the Wolff algorithm. Such updates do 
not change the value of Sj. 

2. The cluster of sites connected to jVoot is constructed in a way analogous to the 
case of 0(n)-spins which was explained in section 3.2. Two sites at positions i 
and j are connected with probability 

pii,j) = max (0, 1 - e-{^[<^r%-C-1-^[*-.<^r]}) , (3.39) 

where 

, -<pf'] - s[<i>r, c'l = ^9ii - j>r>r- (3-40) 

This expression is (up to a factor of two) identical to the spin case, because the 
quadratic contributions cancel. 

3. The phases of the sites j belonging to the cluster are updated according to 

4. If necessary, the whole configuration is shifted by multiples of ±2^ in such a 
way, that the mean value (j) = S^^o^ 'i^j i^ closest to the potential minimum 
corresponding to ji^^^^ = 0. The last step is important to assure detailed balance 
with a finite nmax- The parameter rimax niust be large enough that the re- 
centered configuration is contained in the interval [— nmaxTi", ^^maxTi"]- This value 
can be determined during the thermalization of the system. 
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Fig. 6. Illustration of the cluster algorithm for non-compact phase variables. The figures show 

from top to bottom: (i) Original phase configuration. Possible axis positions are indicated with 
dashed lines, located at the maxima and minima of the cosine potential. The randomly chosen 
axis and root site of the cluster are marked with the black solid line and black dot respectively, 
(ii) Cluster of connected sites. The sites which could potentially connect to the root site are 
indicated with tic marks, (iii) New phase configuration obtained by flipping the cluster around 
the axis. 
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Fig. 7. Integrated autocorrelation times r for {{(f> — 0)^) versus system size A''. The data were 
obtained at the critical point a = 1 using Ej/Ec = 1. 

There are alternatives to fixing the interval of symmetry points and re-centering 
the new configurations, which obviously satisfy detailed balance. One such possibility 
would be to choose the axis among the nmax symmetry points above or below ■ 
Another idea, proposed in Ref. 32), is to choose the axis closest to the phase variable 
of a randomly chosen site. This has the advantage that the axis cuts through the 
phase configuration with high probability, which means that few freezing clusters 
(cluster which span the entire system) in the delocalized phase are produced. 

Using the ideas of Luijten and Blote^'') discussed in section 3.2.3, the cluster 
move can be performed in a time 0(A^log A'^) despite long-range interactions, which 
allows us to compute precise data for large systems (up to A « 10^ if AtEq = 0.25). 
Fig.Elplots the integrated autocorrelation time r for {{(p—^)'^) as a function of system 
size A^. Even though the CPU time for a cluster update is 0(Alog A^) as compared 
to 0(A) for a local update, the gain in sampling efficiency is considerable. 

3.5. Winding number sampling 

In section 2.4 it was shown that the effective charging energy of the single electron 
box can be computed from the winding number distribution of the phase configu- 
rations. If we discretize imaginary time into A time steps At = P/N, the action 
(|2-18jl becomes 



where we have introduced the dimensionless tunneling conductance a = tt-j-d^, 

^ Zn'^ tit ' 

with Rk = ^ (a = ^^^^ for h = 1). Periodic boundary conditions 4)n+i = are 
employed. Since only 27r-periodic functions appear in Eq. H3-42() . (f)j S [0, 27r) can 
be interpreted as the angle which defines the orientation of an XY-spin. Comparing 




(3-42) 
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Eqs. (|2-10j) and (|3-42|) we see that the action of the single electron box is equivalent 
to that of a single dissipative quantum rotor. 
In Fourier space, Eq. H3-42() becomes local, 



N-l 
A;=0 

where V'fc = ^^=q e^~^^e'^'^^ denotes the Fourier transform of e**^ and the Fourier 
transform of the kernel (j 7^ 0) 

Using fast Fourier transformation it is therefore possible to compute the action of 
a configuration in a time O(A^logA^) despite the long-range interactions. This is 
important for the transition matrix Monte Carlo algorithm presented below. 

3.5.1. Path integral Monte Carlo 

The expectation value (w^), and hence the effective charging energy ()2-19() . can 
be computed using path-integral Monte Carlo with cluster updates, as explained in 
section 3.2. This approach works well for small and intermediate values of the tun- 
neling conductance a. Configurations with winding number u; = in are generated 
according to their weight Wn in the partition sum. 



Wn = , ,^ . (3-45) 



However, for large a, the effective charging energy - and therefore also the fraction of 
paths with winding number different from zero ~ decreases as exp(— vr^a) according 
to theoretical predictions. A huge number of configurations with winding number 
zero will be generated for each configuration with non-zero winding number. Due 
to this rapidly deteriorating efficiency, the cluster algorithm can only be used in the 
range of tunneling conductances for which > 10~^Ec- 

3.5.2. Transition matrix Monte Carlo 

In order to compute the effective charging energy at even larger values of the 
tunneling conductance, we developed a new Monte Carlo approach,^) which attempts 
to insert or remove phase slips (kinks) and in doing so measures the relative weights 
of the different winding number sectors. As only paths of winding number and 1 
contribute in the limit of large a, we restrict the discussion to the calculation of the 
relative weight of these two winding number sectors. However, the algorithm can 
also be used to study additional winding number sectors. 

The effective charging energy ()2-19|) can be expressed for large a as 

Eh_ 2vr^ En=on'wn .^ 2vr^ ^1 ... 
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Fig. 8. Illustration of a kink (phase-slip) insertion. The filled dots show the original phase con- 
figuration {(j)n ~ <j>{Tn)\n = 1, . . . , A''} with winding number tj = 0, and the line a phase-slip 
path (/"kink^ periodic- Both the positiou To and the width A of this phase slip path are randomly 
chosen. The open dots show the new configuration with winding number uo = 1, which would 
be obtained by adding (^11°^! periodic to (/>(r). 

The weights wq and wi are proportional to the time, which a Metropohs- Monte Carlo 
simulation would spend in winding sectors and 1, respectively. A "flat histogram", 
which means equal occupation probabilities for both winding sectors, could be ob- 
tained by adding an additional "inverse-weight" factor to the usual Boltzmann term. 
The acceptance probability for a proposed kink update from a configuration with 
winding number n and action S to winding number n' and action S' then becomes 

P(n,S)-.in',s') = min (l, , (3-47) 

whereas a cluster update from n to n' (which takes the change in action into account) 
would be accepted with probability 

=min(l,^ ) . (3-48) 

The "flat histogram" (detailed balance) condition for the a priori unknown relative 
occupation probability wi/wq can then be expressed for n, n' G {0, 1} as 

{Po-.i{wi/wo)) = {Pi^o{wi/wo)) , (3-49) 

where the averages are taken over a random sequence of kink- and cluster updates. 

Our strategy is to determine the transition probabilities (Pq^i) and (Pi^o) for 
different values of wi/wq in order to find the solution of Eq. H3-49() . rather than 
trying to adjust the relative occupation of the winding sectors by a Wang-Landau 
type iterative procedure. This is an approach in the spirit of the transition Matrix 
Monte Carlo method. ^"^^ 

Each winding number sector is treated separately and the cluster updates serve 
essentially to randomize the configurations within that sector, although their con- 
tribution to the probabilities becomes important for smaller a. In a kink 
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update, illustrated in Fig. |S1 we insert a phase slip f/'^lnJ^ ~ arctan ((r — to)/A) 
of random width A at some random position tq, or rather 

,(ro,A) _ arctan((T - ro)/A) 

-^kink, periodic - ±^ arctan(/3/2A) ' ^ ' 

periodically continued outside the interval — ^ < r — tq < ^, which is a slightly 
modified version compatible with the finite size of the system. These are stationary 
paths of the long-range part of the action (|2-18() in the limit /3 — > oo. 

Inserting (j)^^^^'' changes the original configuration (f) with winding number n to cf)' 
with winding number nib 1. The corresponding Boltzmann weights exp(5[(/)] — 5'[0']) 
are computed using Eq. H3-43() and stored during the simulation. Updates which 
produce a configuration with winding number |n| > 1 as well as cluster updates 
which leave the winding number unchanged get the weight 0. A cluster update 
which connects to the other winding number sector gets the weight 1. From these 
data, one can calculate the average transition probabilities {P{i — > j)) as a function 
of wi/wq using Eqs. (|3-47|) and H3-48() . Finally, the relative weight of winding sector 1 
is determined by solving Eq. H3-49() . 

To illustrate this procedure. Fig. |21 shows the probabilities {P{0 1)) (with 
positive slope) and {P{1 — > 0)) (with negative slope) as a function of — ln{wi/wo) 
for several values of a. The intersection points of these curves determine wi/wq. 
The effective charging energies obtained by this new method are consistent with the 
results from cluster Monte Carlo simulations for the values of a which can be treated 
by the latter method. The new approach, however, allows us to simulate the system 
at much higher tunneling conductance. It increases the range of accessible effective 
charging energies by 30 orders of magnitude. 
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§4. Applications 

This section summarizes the main results of our efforts to simulate environmen- 
tally coupled quantum systems, such as dissipative spin chains or resistively shunted 
Josephson junctions. 

In Ref. 1), we studied the phase transition in a chain of Ising spins, using the 
cluster algorithm detailed in sections 3.2 and 3.3. Since the two spin states could 
represent for example the two flux states in a SQUID, such a chain can be thought 
of as a model of a qubit-register. Our careful investigation showed, that no locally 
critical phase exists even for large dissipation strength, and that the second order 
phase transition from the disordered to the ordered state is controlled by a single 
fixed point. In particular, the dynamical critical exponent takes the value z ~ 2. This 
finding, as well as the values of the critical exponents r] and ly, were in remarkably 
good agreement with the predictions from a dissipative ^^-field theory.^)' 

A similar agreement with analytical results was found for the chain of dissipa- 
tive XY-spins. This chain was studied in Ref. 2) as a model for a nano-wire formed 
out of superconducting grains. Besides the bulk phase transition, we studied the 
properties of finite, open chains. If normal leads are attached, the (infinite) wire 
becomes insulating, while superconducting leads turn the whole device supercon- 
ducting. For mixed leads (normal and superconducting), the wire is metallic and its 
zero-bias conductance universal, in the sense that it does not matter how strongly 
the superconducting lead fixes the phase of the neighboring grain. 

In Ref. 4), we have tested the new algorithm described in section 3.5.2 and 
found that it allows to calculate the effective charging energy up to large values 
of the tunneling conductance a. The leading exponential suppression, Eq/Ec ~ 
exp(— vr^a), can be traced over more then 34 orders of magnitude - compared to 2 
(8) orders of magnitude for path-integral Monte Carlo simulations with local (cluster) 
updates. In the zero-temperature regime, we measured a pre-exponential factor ~ 
a^, which is not consistent with any of the numerous theoretical predictions.^-^^'^^^"^^^ 

With the powerful cluster algorithm described in section 3.4 and Ref. 5), we 
were able to verify that the superconductor-to-metal transition in a single resistively 
shunted Josephson junction occurs for Rg = Rq, independent of the Josephson cou- 
pling Ej. Furthermore, the temperature dependence of the resistance in the T = 
superconducting phase is proportional to T'^^^q/^=~^\ as predicted in Ref. 41). Re- 
markably, on the phase transition line, we found continuously varying correlation ex- 
ponents. The exponent r]{q), defined through {exp{iq{4>T — (po))) ~ 

-Ml) (where q 

is some non-integral real number), turned out to decrease exponentially with increas- 
ing Josephson energy. Hence, it appears that the superconductor-to-metal transition 
is in fact a line of fixed points with continuously varying exponents. 

In the two-junction model discussed in Ref. 40), such a line of fixed points 
was predicted to control the superconductor-to-metal transition. Our numerical 
calculations in Ref. 6), which utilize adapted versions of the cluster moves discussed 
in section 3.4, showed, that the critical properties in this interacting system are, 
within error-bars, the same as in the single junction and a mean field theory (which 
assumes non-interacting junctions) was even capable of accurately predicting the 
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phase boundary for intermediate Josephson coupling. These examples illustrate how 
the ability to simulate single junctions and extended systems with high precision has 
led to new insights and will be important in future investigations. 
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